Neural networks determination of material elastic constants and structures in nematic complex fluids

Supervised machine learning and artificial neural network approaches can allow for the determination of selected material parameters or structures from a measurable signal without knowing the exact mathematical relationship between them. Here, we demonstrate that material nematic elastic constants and the initial structural material configuration can be found using sequential neural networks applied to the transmmited time-dependent light intensity through the nematic liquid crystal (NLC) sample under crossed polarizers. Specifically, we simulate multiple times the relaxation of the NLC from a random (qeunched) initial state to the equilibirum for random values of elastic constants and, simultaneously, the transmittance of the sample for monochromatic polarized light. The obtained time-dependent light transmittances and the corresponding elastic constants form a training data set on which the neural network is trained, which allows for the determination of the elastic constants, as well as the initial state of the director. Finally, we demonstrate that the neural network trained on numerically generated examples can also be used to determine elastic constants from experimentally measured data, finding good agreement between experiments and neural network predictions.


Results
The developed neural networks-based method for determining elastic constants is based on the combined modeling of (i) liquid crystal effective dynamics, (ii) light transmission, and (iii) supervised machine learning, as then applicable both to experimental or modeling data. The method starts by calculating a large number of timedependent light beam transmittance functions I(t) during the relaxation of the NLC sample that correspond to different elastic constants, which we then use to train the neural network that recognizes elastic constants from a simulated measurement. Later, the well-trained neural network is used to predict elastic constants also from signals measured from real samples in the laboratory. A general overview of the method, which could, in principle, be used for any other experimentally relevant setup and determining other material parameters, is shown in Fig. 1. Simulation of LC profiles is explained in Methods.
Neural networks based method for determining elastic constants. The orientational dynamics of a nematic liquid crystal depends on the initial state, rotational viscosity γ 1 , cell thickness D and also on Frank elastic constants. Therefore, if the director field is initially deformed by short pulses of the electric or magnetic field, it will reconfigure back to the equilibrium state after the fields are turned off. Consequently, the intensity of the transmitted light through the liquid crystal cell between the crossed polarizers also varies during the relaxation and its time-dependence therefore indirectly carries information about the elastic constants. So the idea of this method is to use a neural network to identify elastic constants from the time-dependent intensities of transmitted light, regardless of the initial state of the director. However, to perform this task, the neural network must be trained on data, i.e. on many examples of time-dependent intensities I(t) and associated elastic constants. In the setup for the determination of the nematic elastic constants, we assume a specific cell geometry where the liquid crystal is confined in a thin cell of thickness D ∼ 10 µm with strong anchoring which is uniform in x and y directions on both boundaries, so that the director, varies only along the z axis and in time t, n = n(z, t) = (n x (z, t), n y (z, t), n z (z, t)) = (cos φ(z, t) cos θ(z, t), sin φ(z, t) cos θ(z, t), sin θ(z, t)) , where θ ∈ [−π/2, π/2] and φ ∈ [−π, π] are spherical angles. For determining only splay and bend constants, K 11 and K 33 , even more simple 2D director geometry (i.e. director profile variability), n(z, t) = (cos θ(z, t), 0, sin θ(z, t)) , proves sufficient.

Figure 1.
Graphical overview of the method for the identification of unknown parameters of dynamical systems based on numerical simulations and artificial neural networks. In this study, such an approach was used to determine Frank nematic elastic constants from measured time-dependent intensities of transmitted light through liquid crystals. www.nature.com/scientificreports/ Determining K 11 and K 33 . The scheme of the method for determining K 11 and K 33 is shown in Fig. 2. To create a training set, we start with generating a pair of random elastic constants K 11 and K 33 from the uniform distribution in the interval of possible expected values. For example, to predict the elastic constants of 5CB liquid crystal at approx 23 °C, we choose elastic constants from the uniform distribution U(2 pN, 18 pN) as we expect the actual constants somewhere within this interval 29,35,36 . Next, the initial non-equilibrium state of the director n(z, t = 0) = (cos θ(z, t = 0), 0, sin θ(z, t = 0)) is set by generating a random, not necessarily physically meaningful smooth 1D function θ(z, t = 0) for z ∈ [0, D] by quadratic interpolation between a random number of points at different random positions within the interval [0, D]. Having the initial state n(z, t = 0) , a pair of elastic constants and the rotational viscosity γ 1 (e.g. 0.098 Pa s for 5CB at the room temperature 37,38 ), the dynamics n(z, t) is numericallly simulated (see "Methods"). The boundary conditions n(z = 0, t) and n(z = D, t) are determined by the selected anchoring type on each boundary. Different combinations of planar ( θ = 0 ) and homeotropic ( θ = π/2 ) anchoring are tested. From the configuration of the director at each time step, the transmission of light through crossed polarizers and the sample of reconfiguring liquid crystal between them are calculated using the Jones matrix formalism 39,40 . We comment that more advanced light propagation methods could also be used, such as finite difference time domain (FDTD) 41 , but give no qualitative difference for the considered nematic geometries. The values of the ordinary and extraordinary refractive indices, the thickness of the cell and the spectrum of the light source are required to be known precisely for particular liquid crystal material and experimental setting. In this way, the time dependence of the transmittance I(t) is calculated and discretized at, for example, 500 time steps for each pair of corresponding random elastic constants and the initial state Figure 2. Scheme of neural network-based method for determining the elastic constants of splay ( K 11 ) and bend ( K 33 ) deformations. First, the elastic constants are randomly set (1), then the random initial state of the director, n(z, t = 0) , is generated (2), and accordingly, the dynamics of the director n(z, t) is calculated (3) that allows us to simulate the time dependence of the intensity of transmitted light through the sample between crossed polarizers (4). This is repeated 200,000 times to generate training and a validation data set. Once the neural network is well-trained, the experimentally measured intensity I(t) can be used to determine the elastic constants of a real liquid crystal sample. For example, for the determination of elastic constants of 5CB, the cell thickness was D = 10 µ m and the time interval was T = 1.28 s. www.nature.com/scientificreports/ of the director. Repeating this 200,000 times, a set of 200,000 pairs of input vectors X i = [I i (t = 0), I i (t = �t), ..., I i (t = T = 499�t)] and expected (true) output vectors T i = K i 11 , K i 33 , packed as {(X 1 , T 1 ), ..., (X 200000 , T 200000 )} , is obtained, which is later split into a training and a validation data set of lengths 185,000 and 15,000, respectively. The time of the interval T depends on the geometry and rotational viscosity γ 1 , and we set it to such a value that the intensity I(t) saturates and effectively becomes constant. For training, the data is -as usual for neural network training 42 -linearly scaled onto the [0, 1] interval. While the relative intensities I(t) are already limited to this interval by themselves, we transform elastic constants by the "MinMaxScaler" as K i and analogously for K 33 . The training data set is used to train the weights and biases of a dense sequential neural network so that the difference between the predicted output Y and the expected output T iteratively becomes as small as possible. To quantify the difference, the mean absolute error was used as a loss function. The validation set is used to test the neural network's performance for data that was not used for training. Training via Adam optimization algorithm 43 was performed using Tensorflow Keras software 44,45 with batch size 25 and learning rate η = 0.0003 . A neural network with an input layer of 500 neurons and four hidden layers of 500, 400, 250, and 100 neurons and rectified linear unit (ReLU) activation functions and an output layer of two neurons with sigmoid activation functions was used. It turns out that the network's architecture can be substantially varied as long as the total number of parameters (connections between neurons) is large enough (for our study, larger than order-of-magnitude ∼ 10 5 ). As observed, adding more layers or increasing the number of neurons in layers in the described model does not significantly improve its accuracy.
A neural network can be trained to recognize both elastic constants K 11 and K 33 from time-dependent transmittances I(t) through samples with n(z, t) = (cos θ(z, t), 0, sin θ(z, t)) director geometry only if the planar ( θ(z = 0) = 0 ) and the homeotropic ( θ(z = D) = π/2 ) anchoring are used at the opposite cell surfaces. Planarplanar anchoring geometry allows us to determine K 11 only, while in cells with homeotropic anchoring on both surfaces, only K 33 can be determined. This is shown in Fig. 3, where 2D histograms show the number count of points in particular sections of the 2D space of true versus predicted constants for examples from the validation set. For a single type of anchoring on both plates, the equilibrium orientation of molecules is always constant Figure 3. Comparison of predicted and actual elastic constants K 11 and K 33 from the validation data set with n(z, t) = (cos θ(z, t), 0, sin θ(z, t)) director geometry and three different combinations of anchoring. Planarplanar geometry (a) allows for precisely predicting K 11 but not K 33 , using homeotropic-homeotropic geometry (b), K 33 can be determined, whereas in planar-homeotropic geometry (c), one can train a neural network to determine both elastic constants K 11 and K 33 at once. The planar-homeotropic configuration is later also used in experiments. www.nature.com/scientificreports/ along z-axis, either θ(z, t → ∞) = 0 or θ(z, t → ∞) = π/2 , regardless of the elastic constants. Consequently, the transmittance after relaxation, I(t → ∞) , does not depend on them and therefore the information about elastic constants is efffectively embedded in the effective dynamics of the time-dependent I(t). In the planar-planar geometry, the dynamics in the last part of the relaxation, when the deformations are small, and the director is only slightly different from its equilibrium configuration n(z, t → ∞) = (1, 0, 0) , is mainly governed by K 11 . This follows from the fact that in a small-deformations regime, (∂n z /∂z) 2 ≫ (∂n x /∂z) 2 and therefore the Frank-Oseen elastic free energy density ("Methods": Eq. 1) simplifies to f FO ≈ K 11/2 (∂n z /∂z) 2 . This is the reason why only K 11 can be determined in such geometry. In homeotropic-homeotropic geometry, it is just the opposite. The dynamics is governed by K 33 and therefore only K 33 can be determined. However, the equilibrium configuration of the director in a cell with planar-homeotropic anchoring geometry is described by θ(z, t → ∞) = πz/2D if K 11 = K 33 and by a convex or a concave function θ(z, t → ∞) when K 11 > K 33 or K 11 < K 33 , respectively. This results also in the dependence of I(t → ∞) on both elastic constants and for this reason, I(t) carries more information in such geometry and it is possible to determine both elastic constants at the same time. Consequently, such geometry has been chosen to be studied in more detail. As depicted by the red curve in Fig. 6, the mean absolute error of the predicted elastic constants can decrease down to σ K ii ≈ 0.5 pN after 60 epochs of training and batch size 25 using the training set corresponding to the parameters, described in the caption of Fig. 4. Of course, the prediction depends on the accuracy of other parameters (refractive indices, n o , n e , thickness of the cell, D, rotational viscosity, γ 1 ) that should therefore be known as precisely as possible. For example, the cell thickness and the refractive indices explicitly determine the phase difference between ordinary and extraordinary polarization and thereby the magnitude of the light intensity, whereas the rotational viscosity directly scales with the elastic constants and in turn with the I(t) curves in time. In Fig. 4, we show the distributions of the elastic constants determined by the neural network that was trained on the data set corresponding to D = 15.0 µm , n o = 1.523 , n e = 1.744 , γ 1 = 0.20 Pa s and tested for time-dependent transmittances that were simulated in systems with elastic constants K 11 = 11.0 pN , K 33 = 17.0 pN , but slightly modified parameters D, n o , n e , γ 1 . Training a neural network with data corresponding to the inaccurate thickness of the cell D or refracitve indices n o , n e causes the predicted elastic constants to be in the wrong ratio, while the inaccuracy of the rotational viscosity γ 1 causes that both predicted constants are shifted by the same factor.

Scientific Reports
Once the neural network is well-trained from numerically generated data pairs ( X, T ) calculated with the same thickness of the cell D, refractive indices n o , n e and the rotational viscosity γ 1 as used in the experimental www.nature.com/scientificreports/ setup, the trained network can be utilized to predict elastic constants of a real sample from an experimentally measured I(t) as well.
Method validation: prediction of K 11 and K 33 from experimental data. Experimentally, the nematic liquid crystal samples of 5CB and E7 were used in D = 10.0 µm and D = 15.0 µm thick cells, respectively, with strong uniform planar anchoring at the bottom ( z = 0 ) and homeotropic anchoring at the top surface ( z = D ) in both cases. The sample of the NLC between crossed polarizers was illuminated by a lamp with a known spectrum and the transmitted light was measured. First, using transparent electrodes positioned at the boundaries of the cell, the director was randomly deformed within the x-z plane by the pulses of the electric field at times t < 0 and after the voltage was turned off at t = 0 , the director freely relaxed from its randomly deformed initial state to the minimum free energy state. 30 experimental measurements I(t) were done for each material and then normalized to the interval [0, 1] and interpolated to 500 time steps within the same time domain as numerically simulated intensities that were used for training, [0, T = 499 t] . Since the performance of our method does not depend on the initial state of the director -due to the variety of the director initial states used for the training data set -the transmittance I(t) could be measured after the sample was already partially relaxed. This allows us that the measured intensities I(t) can be interpolated to many random intervals [ , + T] , where max(�) = T/10 , to increase the number of inputs X and consequently get more results Y = [K 11 , K 33 ] to achieve better statistics of the predicted values.
The measured time-dependent intensities and distributions of predicted elastic constants for the 5CB and the E7 liquid crystal samples at room temperature (approx. 23 °C) are shown in Fig. 5. In each histogram, there are five different distributions of predictions for K 11 and K 33 made by five differently trained neural networks. The weights and biases of the networks are always initially set to random values 42 Table 1. As illustrated in Fig. 4, the accuracy of the predicted elastic constants depends on the accuracy of different material and geometrical parameters, such as D, n o , n e , and γ 1 . In the used (experimental) setup, the thickness of the cell was measured to ±0.1 µ m accuracy, the refractive indices to ±0.001 accuracy, in agreement with typical values reported in the literature [46][47][48] . Rotational viscosity γ 1 was taken from the literature 37,38,49 .
By increasing the number of training epochs, the mean absolute error of predicted constants for numerical examples from the validation set decreases, as shown by the red curve in Fig. 6, but the predictions from experimental measurements of the transmittance I(t) become increasingly inaccurate, compared to the expected values 51 , K 11 ≈ 11.1 pN , K 33 ≈ 17.1 pN , as shown by blue and green histograms in Fig. 6. Besides that, multiple peaks can emerge in the distributions of elastic constants determined by neural networks. We speculate that the reason is the following. Differences between the measured and the simulated time-dependent intensities are inevitable: There are always inaccuracies in the parameters (D, n o , n e , γ 1 ) that are used to generate the training data, the backflow and the light scattering are neglected in simulations, there is also some noise in the experimental measurements. Therefore, we can only approximately simulate the transmittances I(t). Besides that, it is known that neural networks can become overfitted 42,53 and therefore less general and more sensitive to details and noise during training with the increasing number of training epochs. We observe that as long as the model is general enough (after less than ∼ 10 epochs), the distributions of the predicted constants are broader but centered near the actual constants, while after many ( 20 ) epochs, the distributions of the predictions get narrower but less accurate. For this reason, to achieve optimal predictions, one should not over-train the network when the aim is to predict material parameters from experimentally measured signals that are, in details, always slightly different from the simulated ones. Therefore, to avoid overfitting, we stopped training after 10 or even fewer epochs of training. In Fig. 5 and in Table 1, the predicted elastic constants were determined by models that were trained in 5 epochs.
Determining K 11 , K 22 and K 33 . To extract all three elastic constants K 11 , K 22 and K 33 from I(t), it proves that more complex deformations of the director that include twist are needed to emerge in the measuring cell or geometry. To keep the cell geometry, therefore, the director parametrization n(z, t) = (n x (z, t), n y (z, t), n z (z, t)) = (cos φ(z, t) cos θ(z, t), sin φ(z, t) cos θ(z, t), sin θ(z, t)) with infinitely strong anchoring that gives fixed boundary conditions φ(z = 0, t) , θ(z = 0, t) , φ(z = D, t) , θ(z = D, t) is assumed. The initial state of the director in numerical simulations is then determined by two random functions θ(z, t = 0) and φ(z, t = 0) . The transmittance is again simulated similarly using the Jones formalism, and the training of neural networks is almost the same as for the determination of two constants except for the output vector that is, in this case, 3-dimensional (for the three elastic constants), T i = [K i 11 , K i 22 , K i 33 ]. As it is shown in Fig. 7, we have found out that in hybrid aligned nematic (HAN) cells (panel (a)) and in twisted nematic (TN) cells (panel (b)), determination of K 11 , K 22 and K 33 at once is not possible, while in a tilted twist nematic (TTN) cell with planar ( θ(z = 0) = 0 , φ(z = 0) = 0 ) and tilted ( θ(z = D) = π/3 , φ(z = D) = π/2 ) anchoring (panel (c)), the I(t) apparently carries information about all three constants, that can therefore be determined by a properly trained neural network. The experimental realization of such cells is beyond the scope of this work, but is clearly realizable 54 .

Determination of initial director configurations.
As an alternative use of our method, if all material parameters -including elastic constants -are known, our neural networks methodology can also be used to www.nature.com/scientificreports/ predict the initial director field n(z, t = 0) , from the time-dependent transmittance I(t). Below, we show results for effective 2D director geometry with n(z, t) = (n x (z, t), 0, n z (z, t)) = (cos θ(z, t), 0, sin θ(z, t)).
, that is created in a similar way as in "Neural networks based method for determining elastic constants" section but with fixed elastic constants, we train the neural network with an input layer of size 500 and four fully-connected hidden layers built of 800, 600, 400, 200 neurons with rectified linear unit (ReLU) activation functions and an output layer of size 200 with linear activation function. To train the network, the mean absolute error of the predicted θ(z, t = 0) is minimized.
In Fig. 8, we compare neural network predictions of initial states of the director, described by θ(z, t = 0) , with actual ones for differently complex initial states, where the complexity is quantified by the number of inflection points in profiles. It is shown in panel (b) that this method works well when the actual θ True (z, t = 0) has zero inflection points, while the increasing number of inflection points results in an increasing prediction mean absolute error, but as shown in examples in panel (a) in Fig. 8, the approximate shape can still be determined. Like the method for determining the elastic constants, this one is also sensitive to the inaccuracy of the input parameters using which we generate the training set. This is illustrated in panel (c), which shows mean absolute errors of predictions for examples that were generated with slightly different elastic constant K 33 comparing to K 33 = 17 pN used for the generation of training data. Similar discrepancy is expected when other parameters  5 pN, 25 pN) . Panel (g) shows the values of transposed matrices of weights between the last hidden layer and the output layer of five independently trained neural networks that were used to predict elastic constants of the E7 sample. All five neural networks had the same architecture and identical training hyperparameters (learning rate, batch size, number of epochs, etc.), and the same training data were used. This illustrates that there exist, due to a large number of model parameters (weights and biases) and random initialization of them, many neural networks with completely different combinations of weights that still result in very similar outputs of the network. www.nature.com/scientificreports/ (cell thickness, rotational viscosity, refractive indices) used for the creation of training data are different from the actual ones. However, knowing all material parameters, one could use neural networks to determine the initial structure of a nematic liquid crystal from the time-dependent transmittance I(t) measured during the reconfiguration of the NLC to the equilibrium. In experiments, this could allow for the fast automatic determination of the director field that is a result of imposed external fields, including electric, magnetic, or light.

Conclusions
In conclusion, we have presented a method based on machine learning that can be used to determine selected material parameters, specifically the nematic elastic constants. Possible limitations of the method were analyzed, especially, the sensitivity on the values of other -to be known -material parameters. The presented approach of combining numerical simulations and experimental measurements by supervised machine learning and neural networks could also be generalized to determine other parameters of liquid crystals in nematic or other phases, such as Leslie viscosities, birefringence, dielectric anisotropy, anchoring strength, probably using more complex nematic director geometries. There are also no principal limitations for using a conceptually similar approach for determining selected dynamic or static parameters -such as tumbling parameter or degree of order coupling terms -for which established experimental methods are very scarce or do not even exist, including in passive, active, or biological soft matter.

Methods
Nematic and light transmission modeling. The established approach to characterization of nematic orientational order at the mesoscopic scales is by the construction of the total free energy functional 55 . At temperatures below the temperature of the nematic-isotropic phase transition and in the absence of external forces due to external fields or surface anchoring, the leading mechanism that affects the nematic ordering is nematic elasticity with any deformation of the orientational order from the unifrom state increasing the elastic free  www.nature.com/scientificreports/ energy of the system. In the Frank-Oseen formulation, that is based on the director field n , which has n → −n symmetry, the elastic free energy can be written as F FO = f FO dV , where f FO is the Frank-Oseen elastic free energy density and K 11 , K 22 and K 33 are the nematic elastic constants. The K 11 , K 22 , and K 33 terms describe the increase of the free energy due to the splay, twist and bend deformations, respectively. The equilibrium configuration of the director n(r) with minimum full elastic free energy F FO in the absence of external electric or magnetic fields can be found by solving Euler-Lagrange (EL) equations is the molecular field which vanishes in the equilibrium. In principle, free energy (Eq. 1) could also include saddlesplay f 24 = −K 24 (∇ · (n(∇ · n) + n × (∇ × n))) and splay-bend f 13 = K 13 (∇ · (n(∇ · n))) free energy contributions. However, these free energy contributions are relevant only through the boundary conditions 35,56 , and thus in simple geometries can usually be ignorred 57 . Furthermore, especially in nematic geometries, which include topological defects, the tensorial Landau-de Gennes formulation of the free energy is more appropriate to use 40,58 , and actually, the developed method could also be extended to such tensorial modeling of the nematic. We model the relaxation of the nematic from an arbitrary initial configuration to the equilibrium by the simplified relaxational dynamics equations: where γ 1 is the rotational viscosity, which notably ignores material flow and the corresponding backflow coupling [59][60][61][62] . Such simplified dynamics is an approximation, but when compared to experiments, it often proves sufficient in confined systems and simpler geometries, where weak material flow or flow with a simple spatial profile can develop. Due to their anisotropic structure, liquid crystals are optically birefringent, and consequently light which is traveling through them can change its phase, polarization, and direction of propagation. The latter can be Figure 7. Prediction of all three elastic constants from I(t) measured through samples with 3D director with 1D dependence, n(z, t) = (n x (z, t), n y (z, t), n z (z, t)) . Comparison of predicted and actual elastic constants K 11 , K 22 , K 33 from the numerical validation data sets with different anchoring combinations, marked by two-headed arrows in the first column. While in the hybrid-aligned nematic (HAN) cell (a), it is possible to determine K 11 and K 33 , in the twisted nematic (TN) cell (b), only K 11 can be roughly predicted, but using a tilted twist nematic (TTN) cell with planar ( θ(z = 0) = 0 , φ(z = 0) = 0 ) and tilted ( θ(z = D) = π/3 , φ(z = D) = π/2 ) anchoring (c), it is achievable to extract all three constants from the time-dependent transmittance I(t) at once. In this case, the change in polarization of light when passing through a liquid crystal can be described by the basic Jones matrix formalism, where the liquid crystal acts as a phase retarder. If it is placed between crossed polarizers, its director field configuration can affect the phase change and thus the intensity of the transmitted light.

Experimental measurement.
In the experiments, we have used cells with hybrid alignment nematic configuration where the bottom glass was covered with rubbed polyimide (SE-5291, Nissan) to achieve planar alignment and the second glass was covered with DMOAP silane (ABCR GmbH), which ensures perpendicular orientation at the top glass. The thickness of the cells was controlled with Mylar spacers and measured by the standard interferometric method using a spectrophotometer. ITO coated glasses were used that external electric field perpendicular to substrates was applied. The cells were filled with nematic liquid crystals 5CB or E7 using a capillary effect. In the experiment, we have measured the intensity of transmitted light through the LC sample placed between the crossed polarizers using the optical microscope with 20x objective. The easy axis of the planar anchoring at the bottom surface was set at the angle 45° relative to polarizers. The sample was illuminated with two different LEDs with 505 and 590 nm (Thorlabs M505L3 and M590L4). Different starting position of the director profile was achieved by applying 100 ms electric pulse with a random shape using a programmable waveform generator (DG1022Z, Rigol). The transmitted intensity was measured with a photodiode (SM05PD1A, Thorlabs) and amplifier (PDA200C, Thorlabs) in combination with a digital oscilloscope (MS09404A, Agilent) and digital delay generator (DG645, Stanford research systems). The temperature of the samples was kept constant using a home-made heating stage.

Data availability
Sample data supporting this study's findings which are used for training and testing neural networks in the cited Jupyter Notebook tutorial is archived in Zenodo 65 . Complete data sets are available upon reasonable request from the first author J.Z.

Code availability
The method with relevant code is publicly available in an interactive Jupyt er Noteb ook tutorial on Google Colab 63 and archived in Zenodo 64 .